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Abstract 

Empirical modelling often aims for the simplest model consistent with the data. A 
new technique is presented which quantifies the consistency of the model dynamics 
as a function of location in state space. As is well-known, traditional statistics of 
nonlinear models like root-mean-square (RMS) forecast error can prove misleading. 
Testing consistency is shown to overcome some of the deficiencies of RMS error, 
both within the perfect model scenario and when applied to data from several phys- 
ical systems using previously published models. In particular, testing for consistent 
nonlinear dynamics provides insight towards (i) identifying when a delay recon- 
struction fails to be an embedding, (ii) allowing state dependent model selection 
and (iii) optimising local neighbourhood size. It also provides a more relevant (state 
dependent) threshold for identifying false nearest neighbours. 
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1 Introduction 



The construction of mathematical models whose dynamics reflect the obser- 
vations of physical systems is arguably the fundamental task of physics [1]. A 
model's ability to reproduce the observed dynamics is often quantified through 
the distribution of forecast errors. Indeed in traditional time series analysis 
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the "optimal model" is defined as that which minimises the root-mean-square 
(RMS) error [2]. For a nonlinear system, this approach has the undesirable 
property that it may reject the model which generated the data. In this pa- 
per, a new test is introduced which quantifies the level of consistency between 
each iteration of the model and the observations; this test is of particular 
value for nonlinear models where uncertainty in the data due to observational 
noise limits the utility of statistics like RMS error [2,3,4]. Ideally, a model will 
admit trajectories which shadow the entire dataset to within the limits set by 
observational noise [4,5,6]; testing for consistent nonlinear dynamics (CND) 
provides a simple, computationally tractable, necessary condition for shadow- 
ing. The CND approach is applicable to models based on first principles [7,8] 
and data-based empirical models [9,10,11,12,13,14]. The aim is to confirm the 
simplest model consistent with the data, but none simpler. 

As a test of the dynamics, the CND approach quantifies consistency as a 
function of location in model-state space. It is novel in that it (i) identifies re- 
gions where the model is notably imperfect, (ii) can suggest regions where an 
embedding fails (that is, where no deterministic model is consistent), (iii) pro- 
vides a rational for setting a (local) threshold for the method of false nearest 
neighbours [15], and (iv) provides a means of selecting and weighting models 
for multi-model ensemble prediction [4,7]. Models may be either linear or non- 
linear, and the results below are easily generalised to stochastic models [2,16], 
although the discussion here is restricted to the case of deterministic models 
and additive observational noise. While most of the examples below evaluate 
data-based empirical models, the consistent nonlinear dynamics approach is 
applicable to any dynamical model: it evaluates the particular model along 
with background modelling assumptions such as stationarity and the assumed 
form of observational noise. 

Basic goals and methodology are presented in Section 2. Section 3 provides a 
brief and selective review of data-based model building in delay space. Two 
mathematical systems are then used to illustrate the method in section 4 and 
applications to previously published models of three widely studied physical 
systems, each thought to be chaotic, is given in section 5. The implications of 
these results for nonlinear modelling are discussed and interpreted in section 
6 and conclusions are given in section 7. 

[Fig. 1 about here.] 



2 Methodology 

In this section, a constraint on state-dependent forecast error is derived; this 
forms the basic consistency test of CND. A number of simplifying assumptions 
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are made to ease the derivation, often these can be relaxed and more general 
solutions deployed numerically (at the cost of simplicity and computational 
time). 

The most common means for quantifying the quality of a model are based 
on the statistics of prediction errors [2] (for alternatives, see [3,17,18,19]). 
Consider a model F, with observations Sj G IR ms corresponding to system 
states Xj G E™ 1 of the 'true' deterministic system G. Assume m s = m x = m in 
the following 1 . Also assume additive measurement errors t]i (i.e. Sj = Xj + ^J. 
In this case, the one-step prediction error may be decomposed as 



Ep re d — Sj+i — F(Sj) 

= (xi+i + r) i+1 ) - F(xj + rii) 
= [G(xi) - F(xi)] + [r, l+1 + F(xi) - F(x, + Vi ) 
E mo dci ~\~ E no ; se , (1) 

where E mo( jei represents model inadequacy 2 and E no i se represents error due to 
observational uncertainty. Viewing the prediction error E prc d as a sum of these 
two distinct sources of error highlights three facts: (i) prediction errors may 
be due to noise alone, and need not arise from model error, (ii) model error is 
obscured by observational uncertainty, and can only be accurately assessed if 
this uncertainty is taken into account, and (iii) prediction errors may manifest 
extreme fluctuations throughout model-state space due to variations in model 
sensitivity. Note the similarity of this decomposition with that for "model 
drift" used in operational weather forecasting [8,21]. 

[Fig. 2 about here.] 

Prediction error due to additive noise in a perfect model may be expressed as 
a Taylor series: 

E n oisc = rj i+ i - J( x i)*7i - H(x i )»7 i (2) 



where J(x) and H(x) are the Jacobian and Hessian of F(x). The distribution 

1 Note that in general the model-state space will differ from the true state space of 
the system (assuming one exists); the projection from one space to another poses 
several foundational difficulties. This projection operator will be taken to be the 
identity here (for discussion see [4,20] and references therein). 

2 A given model may be inconsistent due to parametric error (inaccurate parameter 
values) or due to structural error (the model class from which the particular model 
equations are drawn does not include a process that might have generated the data, 
given the observational noise[18]). In the second case the model equations themselves 
are inadequate. 
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of E no i se in (2) reflects how variations in the local derivatives of a particular 
model structure affects the local prediction error. When the distribution of 
rj i is known, the distribution of E noise can be determined from (2) and, since 
E pred is known, regions of model-state space with large E mode i can be identified. 
Figures 1 and 2 illustrate this in model-state space where, for simplicity, the 
initial uncertainty is bounded (constrained to be within the circle) and the 
dynamics are locally linear (the circle evolves into an ellipse). The system state 
at initial time Xj and at final time x i+1 are marked by crosses. For model- 
state vectors S; e IR m and observational uncertainty uniformly distributed 
in a sphere of radius e, there exists an associated ball of consistent model- 
states B e (sj) = {$, : ||£ — Sj|| < e}. In practice, only the observed positions Sj 
(pluses) and their associated consistency balls B £ (sj) (circles) are known. The 
consistency of a model's forecast may be tested by evolving the ball B £ (sj) to 
the final time. This set of evolved states, F{B e (s;)}, will resemble an ellipsoid 
if the evolution is locally linear (i.e. when e is sufficiently small and F is 
sufficiently smooth). If F{B e (sj)} and B e (sj+i) intersect, then the prediction is 
internally consistent (Fig. I), otherwise it is deemed inconsistent (Fig. 2). The 
overlap of F{B e (sj)} and B e (s i+1 ) implies the existence of a model trajectory 
which lies inside the ball B e (sj) at the initial time and also inside the ball 
B £ (sj + i) at the final time; that trajectory is consistent with both the model 
dynamics and the observations. 

The most pedestrian version of CND is simply to verify this necessary con- 
dition. The idea is simply to check whether the largest observational error in 
the worst direction is likely to account for the observations. In practice, ob- 
taining F{B e (sj)} by propagating a ball of model-states in IR m may prove to 
be computationally expensive. A first order approximation provides a linear 
consistency check which is easy to implement. The linear approximation of 
(2) is 

Enoise = m+l ~ J(Xi)»7i- (3) 

The first singular value, cr 1 (s i ), of the Jacobian [4,7] describes the largest pos- 
sible (linear) magnification which corresponds to the major axis of F{B e (sj)}. 
Note from Fig. 2 that there is no initial condition consistent under the model 
if the prediction error is greater than the magnitude of the major axis of 
F{B e (sj)} plus the radius of B e (s i+ i). If the noise level is large relative to the 
local curvature of the (model's) solution manifold, then the linear approxima- 
tion may be generalised for locally nonlinear models; for local linear models 
this result is exact. 

The expected value of the observational noise divided by the (local) length 
scale at which quadratic terms become important provides a fundamental 
ratio for each local linear prediction; if this ratio is large one must consider 
local quadratic models. In regions where this ratio is small, and the length 
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scale at which quadratic terms become important is also small relative to 
near neighbour distances, it is advantageous to take smaller neighbourhoods 
thereby improving the local linear model. CND can identify such regions. 

Note that the details of the consistency test depends on the type of state space 
being employed. The state space may be constructed using either multivari- 
ate data or a delay reconstruction of univariate data. In the first case, each 
component of the model-state vector must be forecast and evaluated while in 
the case of the delay reconstruction only one component need be forecast, the 
others generated by a shift operator; this difference will have an impact on 
tests of CND. The type of uncertainty in the measurements may also be either 
bounded (e.g. truncation error) or unbounded (e.g. additive Gaussian noise). 
In the former, the initial observational uncertainty resulting from bounded 
measurement errors lies in a hyper-cube (that is, a box); given Gaussian noise, 
an isopleth of equal probability in the space will be a sphere (or an ellipse, 
depending on the metric adopted). The discussion below is framed in terms 
of quantization error, but this can easily be interpreted in terms of, say. the 
99.9% isopleth of an unbounded noise distribution. 

Many of the applications presented below assume a delay reconstructed state 
space with truncation errors in the observations, implying that 77 is uniformly 
distributed inside a hyper-cube. The diagonal of any m-dimensional hyper- 
cube of side e may be stretched to at most ai(si)y/me and uncertainty in the 
(scalar) verification implies an additional factor of e, so that the consistency 
measure is 

C(Bi,6)= 'f^ 1 ' ! • (4) 
[(T 1 {8 i )y/m+ 1] e 

This definition of linear consistency explicitly states its dependence on both 
the observational uncertainty, quantified by e, and the local stretching due 
to the model structure, expressed through <7i(sj). Since the value of cxi(s;) 
is obtained from the model, C(sj,e) < 1 is a necessary but not sufficient 
condition. If the model's o"i(sj) is too small the inconsistency will be detected, 
but not if ai(sj) is too large. This approach can be generalised to explicitly 
test for overlap between the two regions, and can be extended beyond one- 
step forecasts [3]; nevertheless the simplest case is used here to more clearly 
illustrate the procedure. 

Note that in some regions of model-state space, large prediction errors are 
expected due to large values of Oi(sj). Alternatively, even relatively small 
prediction errors may be inconsistent in a region of state space where <7i(sj) 
is small. Lone inconsistent predictions may be due to outliers [2] in the data 
stream, so significance is afforded by a localised accumulation of inconsistent 
predictions. Such accumulations may be due to failure of the delay space to 
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provide an embedding, erroneous parameter values, or structural error in the 
model. In a nonlinear dynamical system these are intrinsically mixed; and 
they are indistinguishable given only a single model (see [22,23]). Of course, 
regions of model error can easily be located and interpreted when both the 
model and system are known analytically. While this perfect model scenario 
(PMS) is used in section 4 to illustrate the technique, much interest lies in 
physical systems where no perfect model (or model class) is known. Physical 
systems are considered in section 5, where the CND approach is applied to 
three widely studied physical systems using previously published nonlinear 
models. First, the construction of models in delay space is briefly reviewed. 



3 Empirical models 



This section provides a summary of the data-driven modelling techniques 
which are used below. These models provide a mapping F : IR m — > IR, 
Sj i — > §i +T from the present observed state vector Sj to an estimate §i +T 
of the future observed value s i+Tp where r p is the prediction horizon. In the 
discussion in section 2 the model state space was taken to include simultaneous 
observations of the physical variables. In the case of univariate observations, 
Sj will be represented by a m-dimensional delay reconstruction [24,25] 



Sj— (m— l)Td j • • • ? Si—Tdi &i 



(5) 



where is the time delay. As noted above, the use of a delay reconstruction 
allows CND to focus only on the last component of the state vector. 

Local polynomials can be used to provide an approximation of the nonlinear 
dynamics in a neighbourhood B(sj) about a reconstructed state vector Sj. Two 
local polynomial model formulations are employed below: 

m 

F(si) =a + Y, ajSi_ {j _ 1)Td , (6) 
j'=i 

for a local linear model, and 

m mm 
F(Si) = a + a j S i-(j-l)r d + ajm+lSi-(j-l) Td Si-(l-l) Td , (7) 



for a local quadratic model. The model parameters or polynomial coefficients 
a are determined by solving the over-determined system of linear equations 
formed by substituting the k pairs {s K ^, s K ^ +Tp }f =l in (6) or (7) where k(1) 
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gives the indices of the points found in the local neighbourhood B(sj), specified 
either by fixing the number of neighbours k or by choosing a radius r. Both 
local linear (LL) [10,11] and local quadratic (LQ) [13,14] are employed for 
making predictions in this paper. 

Radial Basis Function (RBF) models of low dimensional chaotic systems (see 
[12,13,14] and references thereof) provide a global, empirically based, nonlinear 
system of equations. An RBF model is of the form 

N c 

F(s l ) = Y / ^(\h-^\\) + L(s t ), (8) 
i=i 

where <fi(r) is the radial basis function, £■ are centres, a = {aj}^ contains 
the model parameters and L represents an additional global linear (or higher 
order) polynomial which is often found to improve the approximation [13,14]. 
The two basis functions used in this paper are cubic 0(r) = r 3 and Gaussian 
0(r) = e~ r / 2cr . Following [14] the centres were chosen uniformly in the recon- 
structed state space in regions where the data existed and a was given by the 
average Euclidean distance between centres which were nearest neighbours. 

An attractive feature of the RBF model is that its parameters can be deter- 
mined via a single singular value decomposition [26]. The model parameters 
a = {dj} k j =1 are determined by solving the linear system of equations b = <fra, 
where the design matrix is $jj = 0(| |sj — Cj\ |)> the images are represented as 
b = {sj +Tp }2i and iVj is the size of the learning data set. This is achieved by 
obtaining parameters a which minimise x 2 = 1 1 1> — <fra|| 2 . Both x 2 an d ||a|| 
are minimised by choosing a = <fr^b, where & is the Moore-Penrose pseudo- 
inverse of <fr [26]. Alternatives to using the default weighting (which is, roughly, 
uniform with respect to the invariant measure) are discussed in Section 5.3 
below (see also Smith [14]). In both cases, the solution is linear in the param- 
eters. 



4 Forecasting where a perfect model exists 

In this section models with known inadequacy are used to illustrate that, at 
least within the perfect model scenario, CND correctly diagnoses model in- 
adequacy. Two well studied mathematical dynamical systems are used. When 
forecasting physical systems, of course, there is no reason to believe that a per- 
fect model exists. In such cases, methods like CND must be evaluated based 
on their ability to provide useful information. This is illustrated for a variety 
of physical systems in section 5. In both cases, it is crucial to distinguish the 
model(s) from the system that generated the data (whether the system be a 
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set of equations or a measuring device). 



4-1 Ikeda map 



In the late 1970's, Ikeda [27,28] pointed out that the plane-wave model of a 
bistable ring cavity exhibited period doubling cascades to chaos. Hammel et al. 
[29] then extracted a complex difference equation relating the field amplitude 
at the (i + l) s * cavity pass to that of a round trip earlier. The amplitude x and 
phase y, corresponding to the real and imaginary parts of the field are related 
by 



x i+ i = 1 + fi(xi cos 9 — yi sin 9), 
y i+1 = n(xi sin 9 + ^ cos 9), 

= a-b/(a* + v> + l). (9) 

This map is chaotic for parameter values a = 0.4, b = 6.0 and /j, = 0.9 [29]. 

A family of imperfect models of this system was developed by Judd and Smith 
[18] by considering the Taylor expansions of the sinusoidal functions in (9). 
The translation 9 = —n + uo is employed [18] to yield an expansion about — 7r 
which is near the middle of the range of 9. Neglecting terms of sixth-order and 
higher yields the model equations: 



cos 9 = cos(-n + u) = -1 + uj l - u; /24 



3 la , .5 



(10) 



sin 9 = sin(-7r + u) = -u + uj a /Q - uf j 120. 



Adopting the Ikeda map (9) as the system, time series were generated for 
both the x and y variables with additive noise. The measurement errors were 
uniformly distribution on [— e, e] with e = 0.01. Both a third-order truncated 
model and a fifth-order truncated model were used to make one step ahead 
predictions of the noisy bivariate time series. These predictions were coloured 
grey when consistent and black when inconsistent (Fig. 3). The circle x 2 +y 2 = 
— 1, corresponding to 9 — —n and uj = (where the truncated models are 
exact) is also shown. As expected the models provide consistent predictions in 
the vicinity of a band around the circle and the width of this band increases 
with the addition of higher order terms in (11). 

[Fig. 3 about here.] 

It is often the case that no sufficiently relevant model equations are available 
and that an empirical, data-driven approach is required. There are two com- 
mon strategies in this case: to construct a local model as each prediction is 
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required or to construct a single, global model. For clarity, only local linear 
models and local quadratic models are employed below, while the empirical 
global models used are based on radial basis functions; both approaches were 
introduced briefly in section 3. 

CND considers the local neighbourhood around the base point of the predic- 
tion. In this example, a linear model is fit to describe the evolution of historical 
base points within this neighbourhood towards their future images. This lin- 
ear model is then applied to the base point to yield a prediction of its future 
image. The optimal size of the local neighbourhood B(x, r), at a point x with 
radius r, depends on (i) the data density in B(x, r), (ii) the magnitude of the 
neglected nonlinear terms in B(x, r) and (iii) the measurement uncertainty of 
points in B(x, r) (for a discussion, see [30,31]). Local linear predictions of the 
noisy measurements of both x and y variables of the Ikeda map using r = 0.15 
and r = 0.2 are illustrated in Fig. 4. The fraction of inconsistent predictions 
were fi nc = 0.025 and fi nc = 0.27 for r = 0.15 and r = 0.2 respectively. The 
predictor with the smaller radius of r = 0.15 produces fewer inconsistent pre- 
dictions because the local linear model is better able to capture the dynamics 
in this neighbourhood. This makes sense, in that for larger values of r there 
are more locations where the local quadratic terms play a significant role and 
thus the linear model (which is blind to these effects) fails to yield consistent 
predictions. 

[Fig. 4 about here.] 

The consistency analysis can also be used to explore whether or not there are 
self-intersections within a reconstructed state space. For the Ikeda map, there 
are many self intersections using a reconstruction dimension of m = 2 (Fig 
5a). These self-intersections imply that nearby points, often called false neigh- 
bours, (Fig. 5b) will have images in different parts of the attractor (Fig 5c). 
The singular value decomposition of the local linear map provides a means 
of calculating the expected region (an ellipse under the linear approximation) 
where the images should fall. In this case of self-intersections, the consistency 
analysis correctly shows that the predictions should lie between the two ex- 
tremes of the images, coloured red and blue (Fig. 5c) and the resulting incon- 
sistent predictions reflect a failure of this reconstruction due to m not being 
sufficiently large. 

Schroer et al. [22] have shown that the fraction of the attractor which yields 
poor predictions due to self-intersections scales as r m ~ Dl where r is the size of 
the neighbourhood used for constructing the local constant predictor and D 1 
is the information dimension. A noise-free time series of the x-coordinate of 
the Ikeda map (9) significant to four decimal places was generated and used 
to provide a reconstructed state space with m = 2. Local linear predictions 
with a neighbourhood of radius r were analysed for consistency against a 
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noise threshold of e = 5 x 10 -5 using (4). The fraction of inconsistent points 
fi nc decreases with the size of the local neighbourhood radius r according to 
f inc oc r m ~ Dl (Fig. 6), in agreement with the theoretical calculations [22]. Di 
was estimated using the Kaplan- Yorke conjecture [32] 

= D Lyap = k + (11) 



where Di yap is the Lyapunov dimension, Aj are the Lyapunov exponents and 
k is defined such that Y$=i K > and YhH \ < 0- For the Ikeda system with 
parameter values a = 0.4,6 = 6.0 and \i = 0.9, D Lyap = 1.71 (for discussion, 
see McSharry [23]). 



4-2 False nearest neighbours 

The method of false nearest neighbours introduced by Kennel et al. (see [15] 
and references thereof), is commonly used to select a "sufficient" dimension for 
models based on a delay reconstruction; here "sufficient" means large enough 
to resolve the dynamics without self-intersection of likely solutions in the (pro- 
jected) perfect model . Let dij(m) denote the distance between two state vec- 
tors Sj and Sj in a reconstructed state space of dimension to. This method relies 
on identifying the fraction of false near neighbours, that is pairs of points whose 
separations dij(m) and dy(m + l) satisfy A(sj,Sj,m) = dij(m + l)/dij(m) > 7. 
A shortcoming of that approach to the detection of these false near neighbours 
[15] is the arbitrary threshold 7 which determines which pairs are "false". In a 
state space where the dimension m is sufficient to remove all self-intersections 
A(sj, Sj, to) will be dependent on the structure of the dynamics in the neigh- 
bourhood of Sj and Sj and thus is expected to vary throughout the model-state 
space if Sj varies. 

An alternative method for identifying a sufficient dimension to is to make 
local linear predictions using different values of to. If self-intersections exist, 
the fraction of inconsistent points will scale as f inc = r m ~ Dl , otherwise f inc 
should drop to zero for a sufficiently small value of r. In any case 7 should be 
allowed to vary with x; within the local linear model, for example, it can vary 
with sigma\. 

[Fig. 5 about here.] 
[Fig. 6 about here.] 
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4-3 Rulkov circuit equations 



[Fig. 7 about here.] 

A second mathematical system will be used to illustrate how CND can identify 
state space dependent model error, in this case due to incorrect parameter 
values. The system is the set of equations defining Rulkov's circuit [33,34] arc 



x = y, 

y = -x - Sy + z, 
z = l[af(x(t)) - z] 



o-y, 



;i2) 



where x(t) is the voltage across the capacitor C, y(t) = y j^i(t) with i(t) 
the current through the inductor L, and z(t) is the voltage across the the 
capacitor C. Time has been scaled by -jj^- The parameters of this system 
have the following dependence on the physical values of the circuit elements 
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VLC 
RC ' 




C_ 



(13) 



The function f(x) is 



0.528, x < -x a 

x{\ X ), X a <C X <C X a i 

-0.528, x > x a 



(14) 



and the control parameter a characterises the gain of the nonlinear amplifier 
around x — 0. The parameters of the circuit correspond to the following values 
for the coefficients in the differential equations (12): 7 = 0.294, a = 1.52, 
5 = 0.534, a = 15.6 and x a = 1.2. 

The model is the same set of equations but with x a = 1.4, thus the model 
is structurally correct, has one parameter in error; all other parameters are 
exactly the same as the system. 

CND is now used to map out regions of inconsistent points in the state space. 
The observations consist of three dimensional time series of x, y and z vari- 
ables, forming a state vector x = (x,y,z), contaminated with additive mea- 
surement noise independently and uniformly distributed on [— e, e]; e = 0.01. 
Numerical integration of (12) using a fourth order Runge-Kutta method [26] 
with a fixed integration step yields a discrete-time map. Given an observed 
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initial condition xo = x(io) and an initial uncertainty e(to) at time to, the 
uncertainty after an arbitrary time r is 



e(t + r) = M(x ,r)e(t ), 



(15) 



where the matrix M(x , r) is the linear propagator defined by 



rto+r 



M(x ,r) = exp 



J 3(x(t))dt 



(16) 



L t 



and J(x(i)) is the Jacobian of the flow given by (12). Taking discrete steps of 
size r = 0.01, the model was then used to provide one step ahead predictions 
of the values of the three-dimensional state vector x = (x, y, z). 

CND analysis of a model with x a — lA can be contrasted with that of a perfect 
model (that is, x a = 1.2) to demonstrate how the approach identifies regions 
of model-state space with systematic errors due to parameter uncertainty. 

The first singular value of the linear propagator M(xo, r) provides a consis- 
tency bound (4) for each prediction (Fig. 7). As expected, the perfect model 
with x a — 1.2 yields no inconsistent points (Fig. 7a and 7c) whereas the model 
with x a — 1.4 gives rise to occasional inconsistent predictions (Fig. 7b and 7d), 
the location of which correctly diagnoses the (known) model imperfection. In 
this way, a CND analysis can be used to identify "synoptic patterns" where 
models tend to fail. This identification does not replace the need for insight 
and physical understanding to determine how to improve the model, rather 
CND merely identifies where the model is most vulnerable; and it can be 
applied to the largest of numerical simulation models. 



5 Forecasting physical systems 

In this section CND is applied to three physical systems; in each case the 
performance of previously published models are contrasted. While the state 
space reconstructions appear adequate to resolve the dynamics, these examples 
demonstrate that other problems arise which are model dependent. The local 
linear is inaccurate when nonlinear terms dominate in the local neighbourhood 
to collect sufficient near neighbours to estimate the model parameters. This 
problem often arises when some regions of the state space are sparsely sampled, 
thereby requiring a large neighbourhood. The RBF model provides a global 
approximation to the dynamics, but is likely to yield a good fit of regions of 
state space with a high data density while neglecting those which are poorly 
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sampled [14]. CND provides a means of identifying and addressing all these 
problems and suggests a method for selectively using different models with 
complementary strengths and weaknesses. 

5. 1 Power output from NH 3 laser 

[Fig. 8 about here.] 

Fig. 8a shows the output power of a NH 3 laser whose dynamics are associ- 
ated with the Lorenz-Haken equations [35]; this dataset, collected by Hiiner 
et al. [36], has been widely studied (see Weigend and Gershenfeld [37], and 
references thereof). The noise level is taken to be below 3 the resolution of 
an (8 bit) analogue to digital converter, thus e = 1. Following Smith [38], a 
Local Linear model (LL la ), was employed to make one-step predictions and 
compared with a Radial Basis Function model (RBFi). Both models use the 
same delay reconstruction parameters (see Table 1). Typical segments of the 
time series and images of inconsistent predictions are plotted for both models 
in Fig. 8. Note that all the inconsistent points occur around the collapses; this 
may result either from the relatively sparse sampling of this region of model- 
state space or a local failure of the embedding [38]. A 2D projection of the 
4D reconstructed model-state space is illustrated in Fig. 8c and 8d, showing 
the origins of the inconsistent predictions for each model. If the reconstruc- 
tion parameters, m and are sufficient to resolve the dynamics, it may be 
possible to improve the consistency by employing a neighbourhood size which 
is suitable for the local linear approximation [30]. 

The problem of low data density around the collapses was addressed by using 
a different LL model (LL^) with a smaller neighbourhood size given by k = 8, 
thereby preventing neighbours with different dynamics from heavily influenc- 
ing the predictions in data sparse regions. This new LL model removes some 
inconsistent predictions at the beginning of the collapses, but adds some new 
inconsistent predictions at the end of the collapses (Fig. 8d and stars in 9b). 

The fact that some predictions are inconsistent for one model, yet consistent 
for another (Fig. 9) suggests that a hybrid predictor using both of these models 
selectively would outperform any one individual model; indeed such a hybrid 
model for this system is presented by Smith [38] based on RMS error criteria. 
One advantage of using a CND criteria instead is that it focuses attention 
on regions of model inadequacy, even when the one-step RMS error may be 
rather small. RBF predictions with large absolute errors tend to correspond 
with small LL lb errors and vice versa (Fig. 9). While the RBF has inconsistent 

3 Except for that due to saturation of the analogue to digital converter, see Smith 
[38]. 
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predictions at the beginning of the collapses, LLij, gives inconsistent predic- 
tions at the end of the collapse, just as the intensity starts to increase again. 
Table 1 provides a description of the models and their prediction results and 
Table 2 gives a comparison of the different models. In particular, while the 
RBF and LLi b models generate 1.13% and 0.48% inconsistent predictions indi- 
vidually, they have only 0.1% in common. Furthermore LLu, provides a better 
complement to RBF than LLi a since there are 0.15% inconsistent predictions 
shared between RBF and LL lfI . 

To establish whether or not an in-sample CND analysis of the learning set 
might contain information on out-of-sample forecasts, a hybrid model was 
constructed. Following Smith [14,30], the RBF centres were used to define 
a Voronoi partition 4 of the state space and the fraction pj of inconsistent 
predictions in each partition j (j = 1, . . . , N c ) was calculated by making in- 
sample predictions of the learning data set. For each base point Sj belonging 
to a partition which had one or more in-sample inconsistent predictions, a 
local linear model (either LLi a or LL15) was used to generate a prediction. 
If the inconsistent regions of state space are clustered in state space and the 
local linear models are complementary to the RBF model, then this hybrid 
should improve the prediction performance. Indeed, both the combined mod- 
els, RBFiLLi a and RBFiLLif, had smaller RMS errors than their constituent 
models and also lower fractions of inconsistent predictions (Table 1). 

[Fig. 9 about here.] 
5.2 Power output from an NMR laser 

[Fig. 10 about here.] 

Figure 10a illustrates a stroboscopic view of the output power of a nuclear 
magnetic resonance (NMR) laser operated at ETH Zurich [39]. The lasing 
particles are Al atoms in a ruby crystal, and the quality factor of the resonant 
structure is modulated periodically. Output power of the laser is reflected in 
the voltage across the antenna (thus allowing for negative values). A local 
linear model (LL 2a ) and a RBF model (RBF 2 ) were constructed following 
Kantz and Schreiber [9]. Details of the models are given in Table 1. The 
model LL2a had a fixed radius neighbourhood of size r = 50 (allowing the 
number of neighbours to vary); this is roughly equal to twice the amplitude 
of the reported measurement error amplitude, e = 25. Inconsistent points are 
shown in Fig. 10a and 10c for the RBF model and in Fig. 10b and lOd for the 
LL 2a model. The inconsistent predictions typically originate from the same 

4 A given Sj belongs to partition j, associated with centre Cj> if dij = mm^^d^ 
where d ik = ||sj - Ck\\- 
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location of model-state space. In this case, reconstructed vectors from regions 
where a large positive observation follows a large negative observation tended 
to be inconsistent. 

This insight from the CND analysis suggests examining the scatter plot of 
the prediction errors (Fig. 11a). This figure shows that errors generated by 
the RBF and LL 2fl models are highly correlated, with both making errors of 
similar sign and magnitude. By zooming into the inconsistent region of the 
attractor, it becomes clear that model inadequacy around the elbow of the 
attractor, shown for the base of the inconsistent points (Fig. lib) and their 
images (Fig. 11c), causes both models to fail. This effect is most noticeable 
for the LL2a model at the elbows of the attractor. In contrast, the ability 
of the LL 2a model to resolve dynamics in a small neighbourhood provides 
consistent predictions in the upper part of the attractor (Fig. lib) whereas 
the RBF model is unable to resolve the dynamics around the two leaves of the 
attractor. 

One method for addressing the failure of the LL 2a model to approximate the 
dynamics at the elbow of the attractor (Fig. lib) is to use a local quadratic 
(LQ) model. A LQ model (LQ 2a ) with neighbourhood r = 50 did not change 
the number of inconsistent predictions and increased the prediction errors (Ta- 
ble 1). By decreasing the radius to r = 25, a new LQ model (LQ 2 &) decreased 
the fraction of inconsistent predictions (from 1.4% to 0.5%) at the elbow at 
the expense of increasing the overall prediction accuracy (Table 1). 

The benefit of using two models together may be seen from the fractions of in- 
consistent predictions in the different scenarios (Table 2). While the RBF and 
LQ 2 6 independently have 1.65% and 0.83% inconsistent predictions respec- 
tively (Table 1), they only have 0.375% inconsistent predictions in common. 
Note that the local models using large neighbourhoods, LL 2a and LQ 2(l , offer 
little compensation to the model inadequacy in the RBF model. 

[Fig. 11 about here.] 

An alternative to changing the structure of the model is simply to alter the 
local neighbourhood size; this can even lead to a meaningful 5 reduction of the 
RMS error. This is illustrated in the next paragraph by taking the inconsistent 
points of model (LL 2a ) and predicting them with a LL model (LL 2 b) that uses 
a fixed number of neighbours given by k = 12. While LL 2 ^ has a slightly higher 
RMS error than LL 2a , it reduces the fraction of inconsistent predictions from 
1.4 to 0.5 (Table 1). 



5 A meaningful reduction in that the dynamics are more accurately reflected as 
opposed to, say, the reduction due to providing the mean of a bi-modal distribution. 
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As noted in the introduction, the appropriateness of each local linear model 
can be determined by examining the expected value of the observational noise 
divided by the (local) length scale at which quadratic (higher order) terms 
become important. CND can suggest where in the model-state space this fun- 
damental ratio is large, and in such locations one must consider local quadratic 
models as above. In regions where this ratio is small, and the length scale at 
which quadratic terms become important is also small, it is advantageous to 
take smaller neighbourhoods thereby improving the local linear model. In this 
case CND successfully identifies such regions. Figure 12a shows the expected 
magnification 6 in a LL 2 6 model plotted against the that of the LL 2a model 
for only those points where the LL2 a model is not consistent. The fact that 
the points generally lie above the diagonal indicates that the LL^f, model is 
(locally) more sensitive to uncertainty than the LL la model which averages 
effects over a larger area. Not only does this increased sensitivity correctly 
reduce inconsistency, it also leads to better (local) predictions, as indicated 
by Fig. 12b, the corresponding scatter diagram of observed prediction error. 
For these points the LL26 model has a lower median error, lower RMS error, 
and predicts 66% of the points more accurately. 

[Fig. 12 about here.] 
5.3 Temperature data from a Fluid Annulus 

Figure 13a shows a time series from a temperature probe in a rotating annu- 
lus of fluid [40,41]. A classic experiment in geophysics which Lorenz [42] cited 
as physical motivation for deterministic aperiodic flow, the annulus consists 
of thermally conducting side walls and insulating boundaries on the top and 
bottom. A temperature difference is maintained between the inner and outer 
side walls, providing an infinite dimensional simulation of the mid-latitude cir- 
culation of the Earth's atmosphere. Following [14], one step ahead predictions 
were made using a local linear model (LL 3 ) and a RBF model (RBF 3 ). See 
Table 1 for details of the models. The data was assumed to have uniformly 
distributed noise of amplitude e = 0.1. The inconsistent points (Fig. 13) reveal 
that the LL model is extremely good (fi nc = 0.68%) whereas the RBF model 
(fine = 2.25%) yields a much larger number of inconsistent predictions. 

One approach for improving the consistency of the RBF model is to force it to 
provide a better fit to the dynamics in regions of state space which were incon- 

6 In a delay reconstruction, the orientation of the forecast error is known a priori. 
The expected magnification, 7, is the stretching in this direction (it need not be 
the first singular value); if the data is noise- free the other directions correspond to 
a rotation. In practice, there is uncertainty in all components; contrasting o\ and 7 
provides information about noise reduction that will be considered elsewhere. 
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sistent in the learning data set. The RBF centres were used to form a Voronoi 
partition of the state space and the fraction pj of inconsistent predictions in 
each partition j (J — 1, ... , N c ) was calculated using in-sample predictions of 
points in the learning data set. Following Smith [14] which re- weighted par- 
titions to get a more uniform error in the model-state space, CND analysis 
provides the information required to construct an RBF model which better 
approximates the more inconsistent partitions (those with large values of pj). 
This is done by providing weights w = such that Wi = l + ap n ^ where 

n(i) is the index of the partition containing Sj, and computing new parameters 
a = (w$) t (wb). This new weighted RBF (hereafter, WRBF) model will have 
a higher RMS error since, by definition, the RMS error is minimised when 
Wi — 1 for all i. While the prediction errors will increase in the partitions with 
small pj, the WRBF can be used to reduce the overall number of inconsistent 
predictions. The value of a controls the balance between increasing the RMS 
error and improving the consistency. Using a = 128, the fraction of inconsis- 
tent predictions was decreased, from 2.25% for the RBF model to 1.07% for 
the WRBF model, whereas the RMS (normalised by the standard deviation 
of the data) increased from 0.47 to 0.89. Many of the large RBF errors were 
improved by the WRBF (Fig. 14a) yielding less inconsistent predictions (Fig. 
14b and 14c). The fraction of inconsistent predictions shared by the RBF and 
WRBF is only 0.586% (Table 2). This small fraction is equal to the fraction 
of shared inconsistent predictions for the RBF and the LL. While it is to be 
expected that the RBF and LL would yield different inconsistent predictions 
because of the disparity between their model structures, this result suggests 
that the weights provide a means of obtaining a complementary version of a 
given RBF model. 

A hybrid model RBFWRBF was constructed using the RBF to make predic- 
tions in all partitions except those where there were one or more inconsistent 
predictions in the learning set. In this case the inconsistent predictions are 
not clustered and so the hybrid gives a medium level of performance, having 
a lower RMS error than the WRBF but a higher RMS error than the RBF. 

[Fig. 13 about here.] 

[Fig. 14 about here.] 



6 Discussion 

A new test for Consistent Nonlinear Dynamics (CND) has been introduced 
and illustrated on a number of examples. The key aim of CND is to quan- 
tify the consistency between a model's dynamics and the observations locally 
throughout the model-state space. Regions of systematic failure indicate states 
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of the system where the model needs improvement, regardless of whether the 
absolute value of the errors in that region are large or small. Similarly, regions 
of large errors which are consistent with the model dynamics and level of ob- 
servational uncertainty should not be counted against the model a priori In 
such regions the observations are consistent with the model forecasts, and the 
model can be accepted. Of course, this acceptance is provisional, the model 
may well fail to remain consistent when the noise level is reduced, or another 
(locally) consistent model may be found to have better error statistics in this 
region and longer shadowing times. In any event, lower one-step forecast er- 
rors per se need not indicate a better model. The errors of the best first guess 
forecast are simply beside the point: the simplest model consistent with the 
data should admit trajectories consistent with the observations. 

Where the model has been consistent, prognostic assessment of the likely ac- 
curacy of predictions can be made in real-time, while in regions of the model 
state space where the model tends to be inconsistent model based estimates of 
predictability (whether analytic, or made with ensemble forecasts [7]) should 
not be relied upon; in such regions historical errors may prove to be of value 
[14]. 

Other insights from the application of CND include: 

(i) The distribution of prediction errors from a nonlinear model is expected to 
show correlation in state space and thus show significant residual predictability 
[30,3]. Equation (2) shows that this is to be expected even with a perfect 
nonlinear model: the local residuals need not be symmetrically distributed 
about the global mean error [3,23]. 

(ii) The original method of false nearest neighbours [15] can be strengthened, 
replacing the arbitrary global threshold used to define "false" by adopting a 
local threshold based on the internal consistency of the dynamics. 

(iii) Noise reduction strategies often use a variational approach [7]; the inter- 
pretation of a variational data assimilation scheme assumes the existence of a 
consistent model trajectory. CND can verify this assumption. 

(iv) The 'optimal linear predictor' is often inconsistent with the observations 
in a systematic manner; CND will detect this. 

(v) All models of physical systems contain structural error; no single best model 
need exist (for discussion, see [5,43,44]). Rather than trying to obtain a single 
optimal model, it may be prove effective to consider ensembles over models 
with different model structures. Ensembles over different initial conditions are 
often utilised to account for observational uncertainty [4,5,7]. Results from 
this paper suggest a complementary method which accounts for structural 
uncertainty by employing ensembles over model structures (see [23,43]). 
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(vi) Systematic changes in the frequency or location of inconsistent points may 
indicate non-stationarity in the underlying process [14]. 

(vii) Given a data stream, selective refinement of points corresponding to 
inconsistent regions is expected to yield a better data base for local models 
than uniform sampling, at least in the noise-free case for local polynomial 
models [30,23,45]. 

[Table 1 about here.] 
[Table 2 about here.] 



7 Conclusion 



A new approach for finding the limitations of dynamical models has been 
proposed and illustrated. By examining the consistency of the local model dy- 
namics with the observed system dynamics, the consistent nonlinear dynamics 
(CND) approach can identify regions of the model-state space where the model 
is systematically inconsistent. Analysis of two theoretical systems and three 
physical systems illustrates that this novel consistency test contains a wealth 
of information about a model's ability to approximate the observed dynamics 
in general, and delay reconstructions in particular. If the delay reconstruction 
does not yield an embedding, then all models will be inconsistent in regions 
where the embedding fails. When an embedding does exist, distinct model 
structures may be preferred in certain regions; CND identifies these explicitly. 
Each analysis yields a direct assessment of individual predictions, taking the 
local properties of the model structure into account, thereby avoiding biases 
arising from model sensitivity. 

Uncovering the failure or success of candidate model structures provides a 
useful discriminator for contrasting different models. The ultimate aim here 
is to find better models, and a better way to define "better" in this context. 
Locally consistent models are expected to allow longer shadowing trajectories; 
a comparison along these lines will be presented elsewhere. 

By adopting CND as a goal, one aims to get the simplest models consistent 
with the data, but none simpler. RMS skill scores can prove a distraction in 
this quest. By addressing the interplay between the nonlinear model structure 
and observational uncertainty in the measurement process, the CND approach 
opens many avenues for both the evaluation and the application of nonlinear 
models. 
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Fig. 1. Consistent prediction: a trajectory of G, with initial and final positions x (x), 
corresponding observed positions s (+), their consistency balls, B e (sj) and B e (sj + i) 
(circles), and images of model-states which are consistent with Sj, F{B e (sj)} (el- 
lipse). The intersection of F{B e (sj)} with B e (sj + i) signifies consistency. 
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Fig. 2. Inconsistent prediction: as Fig. 1. In this case no model trajectory exists 
which is consistent with the observational uncertainty since F{B e (sj)} does not 
intersect B e (sj+i). 
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Fig. 3. Consistent (grey) and inconsistent (black) predictions of the full state space 
Ikeda map using (a) third order truncation and (b) fifth order truncation of (11). 
The black line indicates where the truncated models are exact (8 = —it). 
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Fig. 4. Consistent (grey) and inconsistent (black) predictions of the full state space 
Ikeda map using local linear prediction with neighbourhood radii of (a) 0.15 and 
(b) 0.20. 
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Fig. 5. An illustration of self-intersection in a m = 2 reconstruction of the Ikeda 
map: (a) the reconstructed state space showing observed points within 

a square centred at the base-point (1.05,0.2) being consistent with measurement 
errors of magnitude e = 0.025, (b) a zoom-in showing pre-images of points which 
have large xi + \ (pluses) and points which have small values of Xj+i (dots) and (c) 
the reconstructed state space (xi,Xi + i) with images from (b), corresponding local 
linear predictions (black) and consistency ellipse (grey). 
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Fig. 6. Fraction of inconsistent predictions as a function of the neighbourhood ra- 
dius r used for fitting the local linear model (solid). The dashed line indicates the 
theoretical slope given by m — D\ = 2 — 1.71 = 0.29. 
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Fig. 7. Consistency analysis of the Rulkov circuit. The upper panels illustrate the 
prediction errors for the perfect model (a) imperfect model (b). The grey line is 
the prediction error and the black line is the consistency envelope. Points above 
the envelope are inconsistent. The lower panels are 2D projections of the delay 
reconstruction showing both consistent points (grey dots) and inconsistent points 
(x) for perfect (c) and imperfect model (d). 
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Fig. 8. Predictions of the NH3 laser data. The upper panels illustrate the inconsis- 
tent predictions (x) for models RBFi (a) and LLi a (b). The lower panels are 2D 
projections of the delay reconstruction showing both consistent points (grey dots) 
and inconsistent points (x) for RBFi (c) and LLi a (d). 
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Fig. 9. Consistency analysis of the NH3 laser data showing (a) prediction errors 
for the RBFi and LL^, models, (b) 2D projections of the delay reconstruction 
and (c) the time series. Markers indicate one of four outcomes: (i) neither model 
is inconsistent (dot), (ii) RBFi is inconsistent (square), (iii) LLib is inconsistent 
(hexagon) or (iv) both models are inconsistent (cross). Dashed lines reflect plus and 
minus one standard deviation of the error distribution for each model in (a). 
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Fig. 10. Predictions of the NMR laser data contrasting models RBF2 and LL2 Q - The 
upper panels illustrate the inconsistent predictions (x) for models RBF2 (a) and 
LL2a (b). The lower panels are 2D projections of the delay reconstruction showing 
both consistent points (grey dots) and inconsistent points (x) for RBF2 (c) and 
LL 2a (d). 
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Fig. 11. Consistency analysis of the NMR laser data showing (a) prediction errors for 
the models RBF2 and LL2 a models, (b) return map of Sj versus Sj_ Td and (c) return 
map of Sj+ Tp versus Sj. Markers indicate one of four outcomes: (i) neither model 
is inconsistent (dot), (ii) RBF is inconsistent (square), (iii) LL2 a is inconsistent 
(hexagon) or (iv) both models are inconsistent (cross). 
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Fig. 12. A comparison between models LL2 a {r = 50) and LL2fc (k = 12) for predic- 
tions which are inconsistent under LL2 a : (a) magnification factors 7 and (b) absolute 
prediction errors. Circles indicate points that are also inconsistent in LL2t. 
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Fig. 13. Predictions of the Annulus data contrasting models RBF3 and LL3. The 
upper panels illustrate the inconsistent predictions (x) for models RBF3 (a) and 
LL3 (b). The lower panels are 2D projections of the delay reconstruction showing 
both consistent points (grey dots) and inconsistent points (x) for RBF3 (c) and 
LL 3 (d). 
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Fig. 14. Consistency analysis of the Annulus data contrasting models RBF and 
WRBF: (a) prediction errors for the models RBF3 and WRBF, (b) return map of 
Si versus Si- Td and (c) return map of Si +Tp versus Sj. Markers indicate one of four 
outcomes: (i) neither model is inconsistent (dot), (ii) RBF3 is inconsistent (square), 
(hi) WRBF is inconsistent (star) or (iv) both models are inconsistent (cross). 
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Table 1 

Description of the models used for the three physical systems. Reconstruction pa- 
rameters are dimension m, time delay r^, and dataset sizes used for learning iVj, and 
testing Nt (in 1000s). N c is the number of RBF centres and (p(r) gives the structure 
of the RBF. The LL neighbourhood is defined by the number of neighbours k or 
the radius r. Error is the RMS forecast error normalised by the standard deviation 
of the data and %fi nc is the percentage of inconsistent predictions. 



Key 


System 


m 


Td 


Tp 


A T 


A T 




<j)(r) 


fc 


r 


Frror 


/ojinc 


RBFi 


NH 3 


4 


2 


1 


21 


4 


128 


r 3 






0.0514 


1.13 


LLi a 


NH 3 


4 


2 


1 


21 


4 






32 




0.0381 


0.53 


LLib 


NH 3 


4 


2 


1 


21 


4 






8 




0.0797 


0.48 


RBFiLLia 


NH 3 


4 


2 


1 


21 


4 


128 




32 




0.0376 


0.30 


RBFiLLi 6 


NH 3 


4 


2 


1 


21 


4 


128 


^3 


8 




0.0451 


0.27 


RBF 2 


NMR 


3 


1 


1 


30 


4 


128 


r 3 






0.0207 


1.65 


LL 2a 


NMR 


3 


1 


1 


30 


4 








50 


0.0189 


1.40 


LL 2 b 


NMR 


3 


1 


1 


30 


4 






12 




0.0197 


0.50 


LQ 2a 


NMR 


3 


1 


1 


30 


4 








50 


0.0210 


1.40 


LQ 2 fe 


NMR 


3 


1 


1 


30 


4 








25 


0.0243 


0.83 


RBF 3 


Annulus 


5 


4 


18 


2 


4 


128 


e -r 2 /2a 2 






0.4764 


2.25 


LL 3 


Annulus 


5 


4 


18 


2 


2 






32 




0.4207 


0.68 


WRBF 


Annulus 


5 


4 


18 


2 


2 


128 


e -r 2 /2* 2 






0.8989 


1.07 


RBFWRBF 


Annulus 


5 


4 


18 


2 


2 


128 


e -r 2 /2* 2 






0.5402 


1.27 
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Table 2 

Comparison between two models A and B. Percentage of predictions falling into each 
of four scenarios: (i) Neither A nor B are inconsistent, (ii) only A is inconsistent, 
(hi) only B is inconsistent and (iv) both A and B are inconsistent. 



System 


Model A 


Model B 


(i) 


(ii) 


(hi) 


(iv) 


NH 3 


RBFi 


LLi a 


98.500 


0.975 


0.375 


0.150 


NH 3 


RBFi 


LL lb 


98.500 


1.025 


0.375 


0.100 


NMR 


RBF 2 


LL 2a 


97.675 


0.925 


0.700 


0.700 


NMR 


RBF 2 


LL 2 ft 


98.000 


1.500 


0.375 


0.125 


NMR 


RBF 2 


LQ 2a 


97.600 


1.000 


0.775 


0.625 


NMR 


RBF 2 


LQ 2 6 


97.925 


1.250 


0.450 


0.375 


NMR 


LL 2fl 


LL 2 ft 


98.300 


1.200 


0.300 


0.200 


Annulus 


RBF 3 


LL 3 


97.656 


1.660 


0.098 


0.586 


Annulus 


RBF 3 


WRBF 


97.266 


1.660 


0.489 


0.586 
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